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Many types of bacteria can survive under stress by switching stochastically between two different phenotypes: 
the "normals" who multiply fast, but are vulnerable to stress, and the "persisters" who hardly multiply, but are 
resilient to stress. Previous theoretical studies of such bacterial populations have focused on the fitness: the 
asymptotic rate of unbounded growth of the population. Yet for an isolated population of established (and not 
very large) size, a more relevant measure may be the population extinction risk due to the interplay of adverse 
extrinsic variations and intrinsic noise of birth, death and switching processes. Applying a WKB approximation 
to the pertinent master equation of such a two-population system, we quantify the extinction risk, and find the 
most likely path to extinction under both favorable and adverse conditions. Analytical results are obtained both 
in the biologically relevant regime when the switching is rare compared with the birth and death processes, and 
in the opposite regime of frequent switching. We show that rare switches are most beneficial in reducing the 
extinction risk. 

PACS numbers: 87.18.Tt, 02.50.Ga, 05.40.Ca, 87.23.Kg 



I. INTRODUCTION 

Understanding and quantifying the persistence of bacterial 
populations is of major importance for the efficient treatment 
of diseases. While bacterial persistence was uncovered more 
than 65 years ago [1], conclusive evidence for the underly- 
ing mechanism was only obtained during the last decade from 
laboratory experiments at the single-cell level. It has been 
established that an isogenetic population under identical con- 
ditions can still exhibit two different phenotypes. They are 
clearly distinguished by different rates of cell division: "nor- 
mals" multiply fast, "persisters" do it much slower. For the 
same reason, however, normals are much more susceptible to 
antibiotic treatment, while persisters are highly resilient to the 
antibiotic. An individual bacterium can switch stochastically 
(at a certain rate, often without sensing its environment) be- 
tween the two phenotypes [2] (type-II persistence). 

Systems of two interacting subpopulations, such as normals 
and persisters, have been studied in different contexts in the- 
oretical biology [3-6]. Deterministic models of exponential 
(unbounded) growth were mostly employed, and analysis fo- 
cused on the fitness — the time-averaged net growth rate — of 
the total population, see, e.g., Refs. [7-1 1]. In favorable con- 
ditions, when normal bacteria have a high net growth rate, fre- 
quently switching to persisters is merely a burden, as it de- 
creases the average net growth. If the environment changes 
(deterministically or randomly) between different states, in- 
cluding some which represent adverse conditions for the nor- 
mals, e.g., in the presence of an antibiotic, the same frequent 
switching can become beneficial. In this case, the persis- 
ters uphold a base population size during such a stress phase, 
while normals are heavily decimated. By properly tuning the 
switching rates between different phenotypic states, one can 
optimize the fitness of the total population [8]. For two pheno- 
types and two environments, the average time spent as a cer- 
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tain phenotype should be equal to the average duration of the 
environment in which this phenotype is the fittest one. In more 
complicated models (including phenotype-specific response 
and recovery times upon a change of the environment) one 
still finds that, comparing two (genetic) species, the one with 
switching rates better tuned (in the above sense) outperforms 
the other fitness-wise [10]. 

These are important insights into the role that persisters 
play in a growing population. However, the underlying as- 
sumption of exponential growth is tailored to the description 
of competition among different genotypes trying to establish 
themselves by outgrowing others. Here fitness is instrumental 
to survive in the competition, and a good indicator of a spe- 
cific genotypes' prospects. While such an unbounded growth 
can be realized in vitro, the necessary resources and space in 
vivo are limited. To account for this fact, one should introduce 
models with bounded growth [ 1 2]. In a deterministic (mean- 
field) description, the population will then typically exhibit a 
stable fixed point corresponding to an established population. 
In addition, there will be a fixed point describing an extinct 
population. In reality, population dynamics is a stochastic 
process: an established population is subject to noise com- 
ing from the random character of births and deaths. A rare 
chain of events, where deaths dominate over births, eventually 
drives an isolated established population into the absorbing ex- 
tinction state. Thus for an isolated established population, the 
ultimate goal is survival in the face of intrinsic, and also pos- 
sibly environmental, noise. We suggest, therefore, a paradigm 
shift in the analysis of bacterial phenotype switching by focus- 
ing on the population extinction risk. 

With this motivation, we consider a simple two-population 
system of normals and persisters, possibly in a time-varying 
environment mimicking a phase of catastrophic conditions for 
the population. In a constant environment, a proper measure 
of the extinction risk is the mean time to extinction (MTE) of 
the population, see, e.g., Ref. [13]. We show that a higher 
fraction of persisters exponentially increases the MTE even in 
this setting. With a transient catastrophic phase, a more infor- 
mative measure of extinction risk is the extinction probability 



increase (EPI) because of the catastrophe [14]. Here a higher 
fraction of persisters exponentially reduces the EPI. Therefore, 
when viewed from the perspective of population extinction 
risk, the presence of persisters is always beneficial, providing 
an "insurance policy" against extinction in small communities. 
This should be compared with persisters being a mere burden, 
unless in adverse conditions, when viewed from the perspec- 
tive of fitness. 

The remainder of the paper is organized as follows. In 
Sec. II we set up a simple model that describes the interacting 
populations of normals and persisters. We also introduce, in 
the same section, the pertinent master equation and employ a 
WKB approximation which reduces the master equation to an 
effective Hamiltonian mechanics. We formulate the mechan- 
ical problem that needs to be solved and describe a numeri- 
cal iteration method for dealing with this problem. Sec. Ill 
presents a perturbation theory, based on time-scale separation, 
first for favorable conditions, then including a catastrophic 
phase. There we obtain approximate analytic results for the 
MTE or the EPI, respectively, and for the most probable path 
to extinction, and compare them with our numerical solutions. 
In Sec. IV we contrast the biologically relevant regime of rare 
switching with the regime of frequently-switching bacteria. 
We discuss the main findings in Sec. V. 



II. MODEL AND METHODOLOGY 

We consider a well-mixed two-population system the dy- 
namics of which is described by a continuous-time Markov 
process. The number of "normals" is denoted by n, that of 
"persisters" by m. Normals die at a rate that we set to unity 
throughout, and they multiply at a rate B(l - n/N) per indi- 
vidual. In a stochastic model this corresponds to a finite state 
space, with a maximum number n — N of normal individuals. 
N can be thought of as a number of sites each of which can 
carry at most one individual, or as food resources necessary 
to produce offspring. This dynamics coincides with that of 
infected individuals in the SIS model, with fixed total popu- 
lation size N, unit recovery rate of infected, and an infection 
rate B/N between infected and susceptible individuals [15]. 

We now introduce a persister population whose individuals 
do not multiply or die at all. The populations are coupled by 
normal individuals switching to persisters at a rate a, and per- 
sisters switching to normals at a rate /3 . The ratio of these 
switching rates is denoted F - a/fi. In a mean-field descrip- 
tion, the average numbers of individuals are governed by rate 
equations 



A. Noise and metastability 



n - Bn(\ - n/N) - n - an + jim, 
m — an - j5m. 



(1) 



The rate equations have a trivial fixed point (FP) Fq at n - m - 
0, which describes population extinction, and a nontrivial FP 
Fm at «m = N(l - l/B), Tom = F«m- A viable population 
therefore needs B > 1, when F M is stable, while Fq is a saddle 
point. At the stable FP Fm, the ratio between the population 
sizes of persisters and normals is T. 



Even for large population size, intrinsic noise is crucial, as 
it will ultimately drive the system, residing in the vicinity of 
the deterministically stable FP F M , toward extinction. The 
stochastic system is described by the master equation for the 
dynamics of the probability distribution of population sizes, 

VnJtX 
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Here the Kronecker delta 8„.n prevents transition to a state 
with n = N + 1 . Together with the prescription V n <o im - = 
Vn.nKO and V n >N,m - 0, probability is conserved and limited 
to the stripe («, m) € [0, N] x [0, oo). The extinction probability 
Vo.oit) is described by the equation 
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(3) 



When higher moments are assumed to factorize, the mean- 
field equations (1) are recovered by summation over Eq. (2). 

The stochastic system, as described by Eq. (2), has an ab- 
sorbing extinction state n — = m, corresponding to zero 
eigenvalue and eigenstate 5„,o ;m .o of the transition matrix H. 
All other eigenvalues are negative, hence all other eigenstates 
of the probability distribution decay, and the population goes 
extinct. We assume (and verify a posteriori) that, in con- 
trast to all other nonzero eigenvalues, the eigenvalue with 
smallest nonzero absolute value is exponentially small in the 
system size N. This corresponds to a metastable distribu- 
tion centered around Fm [14, 16-25]. The shape function 
of this distribution, normalized to unity, is called the quasis- 
tationary distribution (QSD); we denote it by n njn . The de- 
cay time of the metastable distribution is T » 1 . An initial 
distribution, describing a viable population, first quickly re- 
laxes to the QSD on a time scale ~ l/(B - 1). Then the 
metastable distribution will "leak" to zero, as described by 
the equations T nj „(t) =* % nm exp(-f jx) [for (n, m) + (0, 0)] 
and Po.oC?) - 1 - exp(-?/r), where T is expected to be expo- 
nentially large in N. Using Eq. (2), the QSD n n . m obeys the 
equation 



>./ T , 



(4) 



and with T exponentially large in N, the right-hand side can be 
approximated by zero. Having found 7t„.,„, one obtains T by 
using Eq. (3): T = l/?Ti,o- One can show (see, e.g., Ref. [19]) 
that T is indeed the mean time to extinction (MTE) when start- 
ing from the QSD. We remind the reader that time is measured 
throughout this paper in units of the death rate coefficient of 
the normal population. 



B. WKB approximation 

When N is sufficiently large, one can approximately solve 
Eq. (4) by a Wentzel-Kramers-Brillouin (WKB) eikonal 
ansatz [16, 26-28] 



n nm = exp [-NS(x,yj\ , 



(5) 



where x - n/N and y - m/N are assumed to be continuous 
variables. Having found S(x,y) in the leading order in l/N, 
the MTE can be calculated up to a pre-exponential factor: 



T=l/?Ei,o*exp[iVS(0,0)], 



(6) 



such that 5(0, 0) plays the role of an entropic barrier against 
extinction. 

Plugging Eq. (5) into Eq. (4) and Taylor-expanding S 
around (x,y) to first order, one obtains, in the leading order 
of l/N, a zero-energy Hamilton- Jacobi equation 



H(x,y,dS/dx,dS/dy) = 0, 



where 

H(x,y, p x , py) - Bx{\ - x) (e A 



1) + x (e" p * -1) 

+ ax (e-"< + ^ -1) + py (e"'- lh -l) 



(7) 



(8) 



is the effective Hamiltonian. The corresponding Hamilton 
equations, 

x = Bx{ 1 - jc) e p > -x e- p - -ax e- p ** p * +/3y e p '" ft , (9a) 



y - axe. 



-Pz+Py 



-fr &-», 



(9b) 



p x = -B(l - 2x) (e"> -1) - (e" A -l) - a (e^* -l) , 
p y = -p(e?*-Py-l), 



(9c) 
(9d) 



describe trajectories of the system in the four-dimensional 
phase space of rescaled population sizes x andy and conjugate 
momenta p x and p y . To determine S(x,y), one can calculate 
the mechanical action accumulated along the proper activation 
trajectory, or instanton, of Hamilton's equations of motion and 
ending in (x, y). 

As the Hamiltonian H does not explicitly depend on time, 
H(x, y, p x , p y ) — E is an integral of motion. In view of Eq. (7), 
the energy E must be zero. One type of motion with E - 
occurs in the invariant plane p x - p y - where Eqs. (9a) 
and (9b) coincide with the (rescaled) rate equations (1). Over- 
all, there are three zero-energy FPs of the Hamiltonian flow: 
(0,0,0,0), [1-1/B,r(l -1/B), 0,0] and (0,0,- lnB,-lnB), 
all of them four-dimensional saddles. The first two originate 
from the mean-field FPs, and we will continue referring to 
them as Fq and Fm, respectively. The third FP, which we 
call F , is the fluctuational extinction point: it appears in a 
broad class of stochastic population models exhibiting extinc- 
tion [17, 20, 21, 29]. Note that all the FPs merge into the 
origin upon approaching the bifurcation point B - 1 . 

As the established population resides around Fm, the instan- 
ton must start at this FP. Now, as we look for 5(0, 0), we need 



to choose between the fixed points Fq and F as the final desti- 
nation. It has been shown that only F a can be reached from the 
region x, y > 0, p x , p y + [21, 24]. The instanton, therefore, 
must be a heteroclinic trajectory which starts at the metastable 
FP Fm at time -oo and enters the extinction FP F a at time +co. 
Finding the MTE, see Eq. (6), demands calculating the action 
S = 5(0, 0) along this heteroclinic trajectory: 



5= dt (pq-#)= dt(-pq-H) 



(10) 



{p x dx + p v dy - H dt ), 



where q = (x, y) and p = (p x ,p y ). In a boundary layer of 
width ~ 1 /N around x — and y — the assumption of large 
population size n, m » 1 breaks down. However, for a suf- 
ficiently large system size N, the contribution of this layer to 
the MTE is subleading in the parameter l/N [25, 30]. 



C. Iterative numerical solution 

The two-degrees-of-freedom Hamiltonian (8) has only one 
independent integral of motion: the energy. It is thus non- 
integrable. Therefore, the instanton can in general be only 
obtained numerically. 

In earlier work, "shooting" algorithms were used to inte- 
grate numerically Hamilton's equations of motion for this pur- 
pose, see, e.g., Refs. [14, 20, 21]. Below (Sec. IIIC) we will 
explain why such an algorithm is not feasible in our case. In- 
stead we adapted an iterative algorithm introduced, in the con- 
text of Hamiltonian field theories, in Refs. [18, 31]. Let sub- 
scripts "M" and "0" label the initial and the final FP, respec- 
tively. We fix a sufficiently long calculation time f max to tra- 
verse the trajectory; it should not be too long in order to avoid 
instabilities in the vicinities of the fixed points. The starting 
iteration numerically integrates Eqs. (9a) and (9b) with the 
momenta fixed at their target values p = p , starting from 
the initial condition q(f = 0) = qM and up to time f max - The 
resulting coordinate curve q(f) is now used to fix the coor- 
dinates in Eqs. (9c) and (9d), leaving a system of equations 
for the momenta, which is integrated backwards in time start- 
ing from p(f = f max ) = p down to t = 0. In each follow- 
ing iteration half-step, momenta (coordinates) are fixed by the 
time-dependent solution obtained in the previous step, and the 
coordinates (momenta) are integrated forward (backward) in 
time, starting from the values at the initial (final) FP and up 
(down) to t — f max (f = 0). We found that this scheme rapidly 
converges to the desired instanton. 

To compute the action, we use the expressions in the first 
line of Eq. (10). The difference between these two versions 
is an easy measure of the numerical accuracy that has been 
reached. 

This algorithm makes it possible to obtain, with little effort, 
the most likely path to extinction and the MTE for a broad 
class of population dynamics models when the target FP has a 
different momentum than the initial FP (as it happens here). 



III. INSTANTON TRAJECTORIES 

A. Close to the bifurcation 

To simplify the algebra, we will restrict ourselves to the 
regime close to the bifurcation point B - 1 where all FPs 
merge, and define the distance to bifurcation 5 = B - 1 «; 
1. As can be checked a posteriori, x, y/T, \p x \, \p y \ ~ 8 or 
smaller. Therefore, exponentials in the Hamiltonian (8) can be 
Taylor-expanded. In addition, we assume that the switching 
from the normals to persisters and back is rare: a, /3 <sc 8 <§: 1. 
Under these conditions, the Hamiltonian (8) becomes 

H(x, y, p x> p y ) - xp x (p x -x + 8)-(ax- j5y)(p x - p y ). (11) 

Here we neglected terms ~ <5 4 , and the term (ax + fiy)(p x - 
p y ) 2 /2 ~ a8 3 . This is consistent if a8 2 » <5 4 , that is, 8 <k 
■\fa. The Hamilton equations read 

x = x(2p x - x + 8) - (ax - Py), (12a) 

y = ax-f3y, (12b) 

Px = ~Px(Px -2x + 8) + a(p x - p y ), (12c) 

Py = -fl(p x -Py), (12d) 

and the zero-energy FPs are (0, 0, 0, 0) (trivial FP, Fo), 
(0,0,-8,-8) (extinction FP, F ), and (<5,F<S,0,0) 
(metastable FP, Fm)- 

It is helpful to rescale all quantities by putting x - 8X, 
y - 8Y, p x - 8Px, p y = 8Py, and t - T / 8. The equations of 
motion become 

^-=X(2P x -X + \)-e(TX-Y), (13a) 

dr 

§-=e(TX-Y), (13b) 

dr 

^ = -Px(Px - 2X + 1) + eT(P x - P Y ), (13c) 

dr 

dP Y 

-L = -e(P x -P Y ), (13d) 

dr 

where e = fi/8. These equations are still canonical with 
Hamiltonian 



happens on the long time scale T ~ 1/e » 1. We formally 
introduce a separate slow time variable T - eT to account 
for this separation of time scales, and consider perturbative 
solutions of the form 



X=X (T) + eX l (T,T') + .. 
Px=Pxo(T) + eP xl (T,T') + 

Y = Y (T') + eY y (T') + ..., 
Py = Pyo(T') + eP Yl (T') + . . 



(16) 



Inserting into the Hamilton equations (13) yields a system of 
partial differential equations in each order of e. Note that, in 
contrast to previous work [20, 21], here the dynamics of fast 
variables (normals) drives the slow variables (persisters). 

In the leading order ~ £°, only two equations remain, 
dXo/dr = Xo(2P xo - X + 1) and dP x0 /dT = -P xo (P x o - 
2Xo + 1). This amounts to the one-dimensional system of 
Ref. [14, 29] close to the bifurcation. The solution must sat- 
isfy the energy constraint hxo - XqPxq(Pxq—Xq+\) - 0, hence 
Pxq = Xq - 1: the projection of the instanton to the X-Px plane 
is a straight line between F M and F a (cf. Fig. 1), and this part 
contributes an action s X q - 1/2 [14, 20, 21]. The solutions for 
Xq and P X q are 



X (T) = 



1 



1+e 7 



Pxo(D 



-1 



+1 



(17) 



where we have arbitrarily fixed the position of the instanton 
along the time axis. 

The slow persister variables appear in the order ~ e 1 , 



^+Yo(T')=TX (T), 

dr' 

dP Y o 

-^ - P Y o(T') = -P X0 (T). 

dT' 



(18) 



On the slow time scale of the left-hand sides, the driving terms 
Xq(T) and P X o(T) change with time only in the narrow region 
\T'\ ~ e <k 1; for earlier and later times they are almost 
constant. Therefore, on the slow time scale they can be de- 
scribed as step functions Xq - 6(-T') and P X q - -6(T'). We 
thus solve dlo/dr' +Yq(T') - F9(-T') by matching solutions 
[with Yq(-oo) = F, Y (+oo) = 0] at V = 0, 



h = H/8 3 = XP x (Px -X + l)- e(FX - Y)(P X - P Y ). (14) 



The action becomes S = 8 2 s, where 



s= (P x dX+P Y dY -hdT). 



(15) 



Yo(T') 



T for V < 0, 

Fe r forr>0. 



(19) 



Similarly, we have dP Y0 /dT' - P Y q(T') = 9(T') [with 
P YQ (-oo) = 0, P Y0 (+°°) - -1], such that 



The rare-switching limit corresponds to e <k 1, and we will 
treat it perturbatively in the following. 



Pyo(T') 



for V < 0, 
for T > 0. 



(20) 



B. Solution in a constant favorable environment 

The leading-order behavior of X and Px, the fast degrees 
of freedom, takes place on the unit time scale T ~ 1 . The 
dynamics of Y and Py, the slow degrees of freedom, however 



The phase trajectory projection to the Y-Py plane forms a rect- 
angle and contributes an area syo — F to the action. To resolve 
the small region \T'\ < e, one would need to include sub- 
leading corrections, which would smoothen the discontinuous 
derivatives of Yq and P Y o at T — 0, round off the trajectory, 
and decrease the action by small terms ~ e. 



The total action in the leading order ~ e° reads 

1 r 
, = - + F. 



(21) 



The MTE of the population becomes, up to a pre-exponent, 



T » exp (n8 2 sq) - exp 



l^+r 



(22) 



In comparison, without persisters the MTE is 
=* exp (N8 2 sxo) - exp (Af<5 2 /2), so the persisters cause 
an exponential increase of the MTE of the population. A 
part of the exponential increase comes simply from an 
increased metastable population size: persisters do not 
compete with normals, so there is no "cost" of increasing 
their population (via F), only a benefit against extinction. 
Therefore, let us compare the MTE (22) with the MTE T ld 
of a single-population system of normals, compensated by 
N — > N(l + T). Both systems then have the same carrying 
capacity K = N8( l+T). The ratio of the MTEs is 



^Id= eXp 



K8T 



2(i +r) 



(23) 



still exponentially large at K8 » 1 and not too small T. No- 
table is the effect of increasing the persister fraction T/(l + T) 
which saturates at large F. Equation (23) does not suggest any 
optimal value of F but the largest possible one; we will dis- 
cuss the relation to other results and the biological context in 
Sec. V. 

Interestingly, persisters contribute an action which does not 
depend on the absolute switching rates a and j3, see Eq. (21). 
It may be surprising that an arbitrarily small but finite pertur- 
bation e > yields an exponential change in the MTE with 
respect to e = 0. This is yet another instance of extinction rate 
fragility [22]. As in other "fragile" population systems, the 
explanation to this counter-intuitive effect comes from a time- 
resolved picture [23]. The effective extinction rate is time- 
dependent. At relatively small times 1 «; T <k 1/e, the ex- 
tinction rate is the same as if the persisters were absent (e = 0). 
At longer times T > 1/e, the extinction rate crosses over to 
its asymptotic value which determines the MTE (22) [23]. 

In deriving Eq. (22), we assumed closeness to the bifurca- 
tion and rare switching, i.e., a, /3 «: 8 <sc 1, or equivalently 
e, eT, and 8 <k 1 ; in particular, implying the upper bound 
F «: 1/e. To obtain the approximate Hamiltonian (11), we 
also had to demand a » S 2 (eF » 8); with hindsight this 
can be lifted: solving the (effectively one-dimensional) fast 
subsystem only employs 8 <K 1, while the ansatz (16) only re- 
lies on time-scale separation e <K 1 . As the small parameters 
8 and e describe unrelated mechanisms, the analytical results 
do not depend (to the given order) on eF » 8. The WKB 
approximation is valid, and the resulting MTE T » 1 is ex- 
ponentially large, if N8 2 (l/2 + F) » 1. For that a minimum 
system size N » S 2 is sufficient, when N~ 1 / 2 «: 5 «: 1 
(QSD width much smaller than the distance between initial 
and target FPs). 

Figure 1 compares the instanton found analytically with the 
numerical solution (see Sec. II C) of Eqs. (13) for a moder- 
ately small e = 0.1. Agreement is reasonably good, and we 
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FIG. 1 . (Color online) Instanton (constant environment, close to bi- 
furcation) for T = 1 and e = 0.1 in several projections. Theory 
prediction (dashed blue) and numerical solution (solid red). 
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FIG. 2. (Color online) Actions of Eq. (15) versus the ratio of switch- 
ing rates T, analytical (21) (solid blue line) and numerical result 
(green marks e = 0.2, red pluses e = 0.1). The error bars were 
obtained by using the original action expression and its integrated- 
by-parts counterpart [see Eq. (10)]. 



checked that it improves, in all projection planes, with decreas- 
ing e. Figure 2 shows that the numerically obtained action ap- 
proaches the theoretical value (21) as e — > 0. The deviation 
also decreases as F goes down, as expected. 



C. Effect of a catastrophe 

What is the effect of a "catastrophe", i.e., temporary adverse 
conditions, on the population extinction risk? For a single 



population, this question was addressed in Ref. [14]. Here we 
find that the presence of a persister subpopulation dramatically 
reduces the extinction probability increase (EPI) caused by the 
same type of catastrophe. 

As in Ref. [14], we will model a catastrophe by setting 
B - during a certain period of time t c . This may mimic the 
effect of a drug that inhibits cell multiplication. The system 
history then differs from the one described in Sec. II A. For 
early times, after relaxation of the system to the QSD, the ex- 
tinction probability still increases with time nearly linearly as 
'Poflit) - 1 - exp(-f/r) =* t/x, where x is the MTE of the sys- 
tem without a catastrophe. At a time to -^ X, when Vo,o - T-qo > 
the catastrophe starts, acting for a duration f c <K X. Compared 
with X, this is a short transient, which may however consider- 
ably increase the extinction probability to the value Vq°q . Af- 
terwards, the system is again described by the (downscaled) 
QSD and continues to decay, while the extinction probability 
increases as Vofl(t) — 1 - (1 - V\ °q) exp[(fo + t c — t)/x]. In 
this setting, the MTE is too crude a measure of the effect of 
the catastrophe: it is dominated by realizations surviving the 
catastrophe, resulting in nearly the unperturbed MTE X. In- 
stead, we measure the influence of the catastrophe by the EPI 
A"Po o = "Pq cf ~~ ^oo ■ Up to a pre-exponential factor it is given 
by ' 



A7>o.( 



-NS C 



(24) 



where S c is the mechanical action accumulated along the in- 
stanton [14], see Eq. (10). While it describes a very different 
quantity, one gets, in the leading order, APo.o from the action 
exactly as one gets 1/t in a constant environment, cf. Eq. (6). 
For Eq. (24) to be valid, in addition to NS C » 1 one has to 
demand that the change of the exponent with respect to the 
constant-environment case is large, N{S — 5 C ) » 1 [14]. 

The instanton itself is obtained analogously to the case of 
time-independent transition rates described in Sec. II B. The 
Hamiltonian now explicitly depends on time: Before and after 
the catastrophe, the system is still described by the Hamilto- 
nian (8). During the catastrophe, the effective Hamiltonian 
becomes 

H c =x (e- p * -1) + ax (*-*•+* -l)+Py (e* - * -l) . (25) 

The instanton trajectory now consists of three connected seg- 
ments: the precatastrophe segment starts at the metastable FP 
Fm and is determined by the Hamiltonian (8); the catastro- 
phe segment is described by Eq. (25); the postcatastrophe seg- 
ment leads to the extinction FP F , again governed by Eq. (8). 
We assume that, after the catastrophe ends, there is still a rel- 
atively large population left (with exponentially long MTE). 
Neither H nor H c depend on time explicitly, therefore on each 
segment, energy is conserved: before and after the catastro- 
phe, H = E = 0, and during the catastrophe H c — E c + 0. 
Furthermore, the phase space points matching the segments 
are fixed by the catastrophe duration t c . In turn, this fixes the 
energy E c . 

Again, we consider the system close to the bifurcation, 5 <K 
1, and assume rare switching, a, /3 <K 5 <k 1, such that before 
and after the catastrophe we have the Hamiltonian (11). We 



expect (and check a posteriori) that x, y/T, \p x \, \p y \ ~ 8 or 
smaller. This leads to 

2 

H c *-xp x + ^-( ax - Py)(p x - p y ), (26) 

where we have kept the same orders as for Eq. (11). 

Rescaling all quantities by 5 as in Sec. Ill A, the Hamilto- 
nian during the catastrophe becomes 
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XP X XP} 
—f- + -f- ~ e(TX - Y)(P X - P Y ), 



with the equations of motion 



^ = X + XP x -e(TX-Y), 
dT 8 

dY 

V=e(TX-Y), 

dP x P x P x 

— --^-- + eT{P x -P Y ), 



dPy 

~dY 



= -e(P x - Py). 



(27) 

(28a) 
(28b) 
(28c) 
(28d) 



The rescaled duration of the catastrophe is denoted T c = 8t c . 
The leading terms in dX /dT and dPx/dT are ~ 1/5 » 1: dur- 
ing the catastrophe the population size decays exponentially 
on the fast time scale. 

To get some insight into the impact of the catastrophe, let us 
consider a numerical solution. To this end, we use the method 
described in Sec. II C, where the equations of motion now 
change from Eq. (13) to Eq. (28) at some time (and back after 
a duration T c ). The result is insensitive to this starting time, 
provided it is sufficiently far from t — and t - f max . Figure 3 
shows several projections of an instanton with and without the 
catastrophe phase, for otherwise identical parameters. In the 
top panels, due to time-scale separation the catastrophe seg- 
ment is nearly horizontal — X and Px rapidly decay, persisters 
are (indirectly) affected much later. The bottom panels show 
that a subpopulation size and its conjugate momentum do not 
change simultaneously. For persisters, first the momentum 
builds up, then the population size drops, as in a constant en- 
vironment (see Fig. 1). For normals, on the other hand, the 
situation has changed; the population size now decays earlier 
than the momentum, this will be explained in Sec. HID. The 
sudden onset and end of the catastrophe is reflected by non- 
smoothness of the instanton (except for the Y-P Y projection). 
"Wiggles" due to nonmonotonic X and P x immediately pre- 
cede or follow the catastrophe segment (we confirmed that 
these are not numerical artifacts). One can see that, after an 
initial decay of the normal subpopulation size, it briefly recov- 
ers, only to be hit all the harder by the catastrophe. Afterwards 
there is a short recovery period caused by influx from persis- 
ters (cf. the Y-Py projection). 

The two-population system with a catastrophe shows a fun- 
damental difference from the single-population case: the in- 
stanton is not only changed during the catastrophe phase, but 
the whole trajectory including pre- and postcatastrophe seg- 
ments is affected. This can be understood via the following 
counting argument. 
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FIG. 3. (Color online) Numerically found instanton for T = 1, £ = 
0.2, and 8 = 0.1, without a catastrophe (dashed blue) and with a 
catastrophe of duration T c = 0.5 (solid red). Green dots mark both 
the start and the end of the catastrophe. 



Imagine we try to match, in a ^/-population system with 
piecewise constant Hamiltonian, the three segments of the in- 
stanton. The 2c/-saddle Fm affords a af-dimensional unstable 
manifold of possible end points of the precatastrophe segment. 
A useful parametrization of this point (where the catastrophe 
segment begins) consists in d - 1 "angles" describing different 
trajectories, and a time-like parameter along the trajectories. 
By matching the catastrophe segment of a given duration, the 
phase space point at its end is then fixed as well. At the other 
end, F a affords a af-dimensional stable manifold of possible 
starting points for the postcatastrophe segment, which can be 
parametrized as above. We thus have d + d parameters at our 
disposal describing the possible points at the end of the catas- 
trophe segment and the start of the postcatastrophe segment. 
Since we have to match them in 2d phase space coordinates, 
this picture does not contradict a unique instanton (although 
there still may be more than one solution). 

For the single-population case d = 1, the phase trajecto- 
ries leaving and entering the fixed points are unique, and thus 
cannot be affected by the catastrophe part in between. In the 
generic case d > 2, however, the pre/postcatastrophe seg- 
ments may differ from the no-catastrophe instanton. In the 
concrete model studied here, these segments have to differ 
simply because of time-scale separation. During the catastro- 
phe, normals are rapidly decimated, whereas the persister dy- 
namics follows much more slowly. In the X-T-projection, the 
catastrophe segment is thus less steep than the slope between 
any two points on the no-catastrophe instanton. It is therefore 
impossible to simply splice the catastrophe segment into the 
latter. 

This explains why "shooting" algorithms are impractical 



for finding the catastrophe-related instanton numerically in a 
multipopulation system. For a single population with catas- 
trophe [ ], such an algorithm can start with a small dis- 
placement from the metastable FP along the no-catastrophe 
instanton, testing different starting points of the catastrophe 
segment — this works as the precatastrophe segment is un- 
changed. Likewise, one can parametrize the zero-energy 
trajectories leaving the initial FP in the two-population sys- 
tem without catastrophe (see Sec. IIIB) by a shooting angle. 
Adding a catastrophe provides an additional freedom (in the 
form of the starting point), and the method is no longer practi- 
cal. 

At the same time, Fig. 3 shows that the instantons without 
and with catastrophe practically coincide (in all projections) 
for an extended part next to both FPs, before eventually depart- 
ing from each other. This means that the system is extremely 
sensitive to minute variations in the angle at which the trajec- 
tory leaves (enters) the initial (final) FP, which only become 
visible closer to the catastrophe segment. We confirmed this 
behavior in tests of the aforementioned shooting algorithm 
(without catastrophe), which, for this reason, already proves 
to be rather tedious. 



D. Analytic theory with catastrophe 

We look for an analytic solution analogously to Sec. IIIB. 
Time-scale separation is still effective: X and Px show fast 
dynamics on the time scale T ~ 1, or even T ~ 8, see Eq. (28). 
They drive the slow Y and Py, which change on a time scale 
T - eT. We denote the catastrophe duration on this scale by 
t; = eT c . 

The leading-order equations ~ e° reduce to the normals- 
only system again, and hx, c - —XPx/8 + XP x /2 governs the 
dynamics during the catastrophe. Since X, Px ~ 1 <sc 1/5, we 
neglect the second term, and arrive at the simple catastrophe 
Hamiltonian hx. c - -XPx/8 used in the single-population 
model [14], The solution is an exponential decay (growth) of 
X (P x ) at a rate 1/5 and for a duration T c . Let X + > X' 
and > P£ > Px denote coordinates and momenta at 
the start and the end of the catastrophe, respectively. Then 
X- = X + exp(-r c /5) and P x = P£ exp(+r c /5). The solution 
forX and P x before and after the catastrophe is the same (up to 
a time shift) as in the constant environment, Sec. Ill B. This is 
no contradiction to the arguments of Sec. Ill C, since the lead- 
ing approximation is effectively one-dimensional. Therefore, 



(l+e r - r < 



forXo >X + , 

(29) 
forXo <X, 

and P X o(T) = X () (T) - 1 for both P X o > PJ = X + - 1 and P x0 < 
Px = X~ - 1. The quantities T < and T > are yet undetermined. 
From the constraints, we get 



X d 



1 



1 + e^/- 5 ' 



P * - i + e ±r c /5 . 



(30) 



and the conserved (X-part) energy during the catastrophe be- 
comes hxo,c - coslT 2 [T c /(25)]/(45). Fixing the time such 



that the catastrophe occurs between T - ±T c /2, we obtain 



Xo(T) 



'I +e T-TAl/S-l/2)\ 
exp(-r/g) 



2cosh[7L/(25)] 

'l +e T+Ul/S-l/2)y l 



forr<-|, 

for 



i < T < +^ 

9 — L — T ? ! 



for | < r. 



(31) 

The momentum is Pxo = -Xo - 1 before and after the catastro- 
phe, and during it decays as 



Pxo(T) = 



- P T / S 



2cosh[r c /(25)] 



= -M-T). 



(32) 



The action found for this Hamiltonian and trajectory is sxo, c = 
[l + exp(r c /5)] [14]. During the catastrophe, the "tra- 
jectory contribution" J PxodXo and J -hxo.cdT cancel each 
other. 

The slow equations of motion (28b) and (28d) are the same 
as in the favorable environment of Sec. IIIB, hence the slow 
leading-order equations (18) (and boundary conditions) are 
unchanged. Again, we only resolve the slow dynamics here. 
The driving terms Xq and Pxo are different now, since a part of 
their movement is replaced by a faster exponential decay (rate 
1/5 » 1) during the catastrophe. Therefore on the slow T 
scale one obtains a step function as an even better approxima- 
tion than in Sec. IIIB. The only difference between Eqs. (17) 
and (31) is that the driving by Xq (Pxo) ends (sets in) at the 
start (end) of the catastrophe T - +T c /2 (instead of T - 0), 
such thatXo = 9(-T^/2 - 7") and P x0 = -0(7" - T£/2). 

Since coordinates and momenta remain separate in 
Eqs. (18), the general piecewise solutions for Yq and Pyo are 
unchanged, but now matched at T - +P//2: 



and 



To(r') 



Pyo(T') 




for 7" < -T c '/2, 

for - r c '/2 < r, 



. e T'-r;/2 
■1 



for 7" < r c '/2, 
forr c '/2< 7". 



(33) 



(34) 



The simple geometric picture that the catastrophe merely time- 
shifts To and Pyo into opposite directions results in a hyperbola 
YqPyo - — Texp(— r c ') on the corresponding segment. 
Persisters contribute an action 



syo,c = / Pyo dPo - 



lYO.c 



d7\ 



(35) 



with the switching Hamiltonian /iy. c = -e(TX - Y)(Px - Py). 
The energy during the catastrophe is evaluated on the slow 
time scale, such that Xq - = Pxo, and 



*70,c 



-E(rXo-Y )(Pxo-PYo) = ere- 



(36) 



The contribution to the action -hyo x T c - -IT c 'exp(-r c ') 
again cancels the phase space area under the catastrophe seg- 
ment, / r exp c Pyo dio- Hence the persister action is syo. c = 
rexp(-r c '), and the total action becomes 



1 



*0,( 



1 + & T C /S 



+ re~ 



(37) 



Reinstating the original time scale t by using 7" 
e8t - j3f we obtain from Eq. (24) 



A'Po.o - exp 



-N8 2 



1 



1 +e'^ 



+ Te 



■pt< 



= eT = 



(38) 



The system without persisters (V - 0) has an EPI =* 
exp (-A^5 2 ixo,c) = exp [-N8 2 / (l + e' c )] . As for favorable 
conditions, we compare with the EPI AVn g of such a single- 



population system of normals, compensated by A^ - 
to have the same carrying capacity K = N8(l + T): 



APq, ( 



exp 



K8T 

i + r 



.-J3'c , 



l 



1 +e'< 



A^(i + r) 



(39) 



The system with persisters has exponentially smaller EPI, to 
which the initial population size K and the persister fraction 
F/(l + T) contribute as to the MTE ratio (23). The paren- 
thesized factor quantifies the fundamental benefit of persis- 
ters and generalizes the numerical value 1/2 in Eq. (23): the 
effect is most pronounced for catastrophes which are long 
on the fast scale of normals, but short on the slow persis- 
ter time scale, t c » 1 » T c ' = j5t c . Then AP ofi /APo d =* 
exp [— K8T/ (1 + T)] , i.e., the ratio is squared with respect to 
the MTE ratio (23) in a constant favorable environment: the 
benefit of persisters is even more apparent in the face of a 
catastrophe. Again the result (39) suggests to choose T as 
large as possible, on which we comment in Sec. V. 

These results are based on 8, e, eF «c 1 (cf. the end of 
Sec. IIIB). For a short catastrophe t c ~ 1 or smaller, the 
WKB result (24) is valid if the reduction N8 2 (sq - «o,c) due 
to the catastrophe is sufficiently large, yielding the condition 
Af » 45~ 2 /f c . A long catastrophe T^ ~ 1 (or larger) strongly 
reduces the action, and the stricter condition is that the remain- 
ing action be large enough. Considering P ~ 1 for simplicity, 
the persister action dominates, leading to N » exp(7 , c ')/(5''r). 

The theory path to extinction is shown in Fig. 4 and com- 
pared with the numerical solution (see Sees. II C, IIIC). For 
a short catastrophe T c = 0.2, persisters are mostly unaffected, 
while the X-Px projection resembles the one-dimensional sys- 
tem [14]. Already for the moderate T c - 1 (not shown), nor- 
mals have gone virtually extinct at the end of the catastro- 
phe, and the population survives mainly due to the remain- 
ing persisters. With a long catastrophe T c - 10, the action 
contributed by persisters is severely decreased as well. Agree- 
ment between analytical and numerical solutions is better than 
in a constant environment. Normals go extinct nearly exclu- 
sively during the catastrophe, which completely determines 
the fast part of the trajectory, rendering the instanton very sim- 
ple. In turn, back-reaction of persisters becomes less impor- 
tant, and replacing the fast driving terms by step functions 
on the slow time scale becomes more accurate. These are 
the main approximations of the zeroth-order theory, hence 
the predictions improve with increasing catastrophe duration. 
We also confirmed that in all projections, the theory becomes 
more accurate with decreasing e. At the same time, the "wig- 
gles" identified in Sec. IIIC become less pronounced. Both 
tendencies go hand in hand, as both are based on reducing 
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FIG. 4. (Color online) Instanton for T = 1 , e = 0. 1 , and 5 = 0.1, with 
a catastrophe of duration T c = 0.2 (top), 10 (bottom), respectively. 
Prediction by theory (dashed blue) and numerical solution (solid red). 
Green dots mark both the start and the end of the catastrophe. 



back-reaction. In Fig. 5, we compare the action (37) with nu- 
merical results. Even for moderately rare switching (e = 0.1), 
the analytical prediction is reasonably accurate, and improv- 
ing with increasing catastrophe duration. 

We summarize the effect of the catastrophe in the leading or- 
der of rare switching. Independent of its duration, the strength 
of the catastrophe is set by the (normalized) death rate of nor- 
mals. Normals decay exponentially on the very fast scale 
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£ = 0.1, and S = 0.1. Error bars span the results obtained using the 
original action expression and its integrated-by-parts counterpart. 



t ~ 1, responsible for the major part of phase space motion 
(unless t c <K 1). Persisters are affected indirectly via switching 
between the two populations. For a short catastrophe, !T C ' <K 1, 
the effect on Yq and Pyq is negligible: switching hardly occurs 
during T c ', and the slow dynamics cannot resolve the differ- 
ence in driving. Therefore only the normal action is reduced, 
and persisters are perfectly buffered against the catastrophe. 
Note that the time f c <K l/(8e) can be much longer than the 
typical lifetime of an individual normal ~ 1. If the catastrophe 
is long enough to be seen on the slow scale, r c ' ~ 1 or larger, 
switching has an effect. While persisters still cannot resolve 
the accelerated extinction of normals, they trace the delay be- 
tween Xq and Pxo in the instanton. On the slow switching time 
scale it appears far shorter, forming a buffer that mitigates the 
catastrophe. The structure of the EPI (38) is thus based on 
the separation between the time scale of the catastrophe effect 
(strength), and the far longer time scale of persister dynamics. 
The catastrophe affects both populations, but acting on nor- 
mals, its duration is measured on the very fast scale t ~ 1 of 
the death rate [action scales ~ exp(-f c )]; acting on persisters 
it is measured on the slow scale T ~ 1 of switching back to 
normals [~ exp(-r c ')]. The crossover shows prominently in 
Fig. 5. 



IV. WHY ARE SWITCHING RATES SMALL IN NATURE? 

So far, we have considered small switching rates between 
the normal and persister states, e <k 1. The correspond- 
ing time-scale separation was the basis of our qualitative ex- 
planation and analytical treatment of the system's dynamics. 
We also numerically examined what happens at e ~ 1 or 
larger. We still consider the system described by the Hamil- 
tonians (11) and (26), respectively, as motivated at the end of 
Sec. IIIB. The instanton and the associated action are again 
obtained as detailed in Sees. II C, IIIC. 

We found that both with and without catastrophe, instan- 
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ton trajectories are qualitatively similar to the e «: 1 case 
even when e = 1 . Further increasing e "locks" persisters ever 
stronger to the dynamics of normals, see Eqs. (13) and (28). 
For very large e, Py =* P x and Y =* TX with only small devia- 
tions. Moreover, persisters still increase the action compared 
with a normals-only system of the same carrying capacity. We 
examined the action as a function of varying switching rate e 
and catastrophe duration t c (N, 8 and F, and hence the car- 
rying capacity K, being fixed). As expected, for given e the 
action decreases with increasing catastrophe duration t c , and 
this decrease becomes stronger for larger switching rate e: the 
more frequent the switching is, the less insurance against ex- 
tinction persisters provide. For given t c , the action decreases 
with increasing switching rate e, and this decrease becomes 
stronger for longer catastrophe duration: persisters are espe- 
cially beneficial in the face of a catastrophe. 

For very frequent switching, there is a new time-scale sep- 
aration which permits an analytical treatment. Consider the 
case 8 «c 1 «: a, j3, (e, er » 1/8), such that switching is 
frequent compared with the normal dynamics even during the 
catastrophe. In both the favorable [see Eqs. (13)] and catas- 
trophic [see Eqs. (28)] environments, we have 



Y = FX- 



1 AY 
edT' 



Py — Px 



IdPy 
e dT 



(40) 



For large e the second term is a small correction, and we ob- 
tain 



FdX 

Y =YX + 

edT 



P Y = Px + ~^r + .-.. (41) 
e AT 



Inserting into the normal equations of motion yields, in the 
leading order in 1/g, the normals-only equations, but with a 
rescaled time T - T /(l + T): 



^=X(2P X -X + 1), fS 

dr dr 

in a favorable and 



-P X (P X -2X+1) (42) 



dX 


X 


dP x 


Px 


=z 


—F +XP X , 






dy 


8 


dr 


8 



—£■ = -§- - -£ (43) 



in a catastrophic environment. As in Sec. HID, from this we 
getX as of Eq. (31), only with the substitutions 7( C ) — > 7( C )/(1 + 
F), and likewise for the momentum P x and the energy h Xx . Y 
and Py are given by Eqs. (41). 

Calculating the action along this instanton, first note that 
the switching term in the Hamiltonian is fty(, C ) ~ 1/e at all 
times. Second, the corrections in Eq. (41) do not contribute to 
the leading-order action, which becomes 



s c * (P x dX + P Y dY) + (P x dX+PydY- h Xfi dT) 

pre /post cat. 

f f (44) 

- (1 + H / Px dX + (1 + O / P x dX - h XiC T c . 

pre /post cat. 

The second and third terms cancel; factoring out (1 + F), both 
contributions are the same as in the rare-switching case, only 



with the above rescaling applied to all times. The first integral 
is also known from the rare-switching case, where it coincided 
with the total action contributed by normals. Applying the 
time rescaling, the action thus becomes 



1 +r 



1 + e r c /[S(i+n] ' 



(45) 



We confirmed (for £ = 100 and various values of T and t c ) 
that this agrees excellently with the action found numerically 
as described at the beginning of this section. This result is 
easily interpreted; very frequent switching effectively "mixes" 
the two subpopulations, as they rapidly switch back and forth. 
Compared with a normals-only system, the factor 1 + F in 
the numerator describes the increased size of the combined 
population. A more subtle effect is the reduction, by the same 
factor 1 + F, of the effective duration of the catastrophe. This 
reduction accounts for the lag still gained by switching to the 
persister state. 

In a favorable environment (f c = 0), persisters switching fre- 
quently do not provide any benefit compared with a normals- 
only system of the same carrying capacity. With a catastrophic 
phase, however, we obtain 



A7>o.o 



An'o 



7= ex P 



-K8 



1 



1 



1 + e'./d+O l + e 'c 



(46) 



This is still a substantial benefit, although much less (for F not 
too large) than that for rarely switching persisters, see Eq. (39). 
Note that here F only appears in the effective catastrophe du- 
ration, not as the persister fraction. 

Persisters are thus most valuable when stochastic switching 
is relatively rare. The fact that rare switching dominates in 
nature can be attributed to an evolutionary process. 



V. DISCUSSION AND CONCLUSIONS 

We have used a simple two-population model of normals 
and persisters to show that (and how) persisters exponentially 
decrease the extinction risk of an established bacterial popula- 
tion. We have compared the two-population system of nor- 
mals and persisters to a normals-only system starting from 
the same total population size. Already in a constant envi- 
ronment favorable for normals, it is beneficial to switch to the 
persister state: persisters contribute to the MTE exponentially 
more than normals since their extinction is delayed by first 
switching back. When the population is under stress — that 
we model as a catastrophe — the same buffering is effective, 
rendering persisters far less prone to extinction, so that they 
exponentially reduce the EPI due to the catastrophe. For catas- 
trophes which are long compared with the lifetime of normals 
but short compared with the much longer switching time scale 
(from persisters to normals), the reduction factor is squared 
with respect to the MTE increase in a constant environment: 
persisters are even more valuable for the population if it faces 
a catastrophe. 

In exponential-growth models which focus on fitness, per- 
sisters are only a burden in a constant favorable environment. 



11 



To explain their existence with an overall benefit one needs to 
invoke temporary adverse conditions. In contrast, we have 
shown that persisters are always beneficial as an insurance 
against the extinction of an established population, as mea- 
sured by the increased MTE, or by the reduced EPI during 
a catastrophe, respectively. We have also shown that to pro- 
vide the optimal benefit, switching to and from the persister 
state has to be rare compared with all other processes. In this 
sense, the extinction risk perspective presented here explains, 
in a natural and robust way, the existence of persister pheno- 
types in bacteria as well as the small switching rates from the 
normals to persisters and back. 

Our main analytical results (23) and (39) advocate that 
switching back from persisters to normals should be rare com- 
pared with switching to the persister state, leading to the 
largest possible fraction of persisters in the metastable state 
(within the range where our theory applies). For a bacterial 
population optimized solely against extinction from the estab- 
lished state, this would be an intuitive strategy even in favor- 
able conditions. During the growth stage, on the other hand, 
the population needs optimal fitness to establish itself. These 
two complementary strategies, optimizing two different quan- 
tities, are not incompatible: the extinction risk perspective ex- 
plains the mere existence of persisters, already without invok- 
ing environmental variations. The switching rates themselves 



(fixing the metastable persister fraction) may be tuned by evo- 
lution to optimize the growth stage in a variable environment. 
Our simple model neglects many features that can be bio- 
logically relevant. For example, in reality persisters have re- 
duced but nonzero birth and death rates. Such a more realistic 
system still features time-scale separation; persisters are now 
directly affected by a catastrophe (e.g., inhibiting their births), 
but again on a much slower scale than normals. Therefore we 
expect a qualitatively similar behavior. Future work can at- 
tempt to account for the cost of switching to persisters, for ex- 
ample, via competition between persisters and normals. There 
are also many alternatives for the detailed dynamics during the 
catastrophe. For many of them, the WKB approximation to 
the master equation provides a viable theoretical framework 
for determining the long-time behavior of bacterial popula- 
tions. 
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